UNCLASSIFIED 


_ AD  NUMBER _ 

AD120083 

LIMITATION  CHANGES 
TO: 

Approved  for  public  release;  distribution  is 
unlimited. 


FROM: 

Distribution  authorized  to  U.S.  Gov't,  agencies 
and  their  contractors ; 

Administrative/Operational  Use;  20  JAN  1956. 
Other  requests  shall  be  referred  to  Office  of 
Naval  Research,  Washington,  DC  20360. 


_ AUTHORITY 

onr  ltr  9  nov  1977 


THIS  PAGE  IS  UNCLASSIFIED 


DISTRIBUTION  STATEMENT  A 


PM  PUBLIC  RELEASE; 


DISTRIBUTION  UNLIMITED 


HAS  BEEN  DELIMITED 


AND  CLCANCD  PM  PUBLIC  RELEASE 
UNDER  OOP  DIRECTIVE  5200.20  Af 
RE  STRICT  I  MS  ARE  IMPOSED  UPM 
ITS  USE  AND  DISCLOSURE. 


f 


Reproduced  by 

DOCUMENT  SERVICE  CENTER 

KNOTT  BUILDING,  DAYTON,  2,  OHIO 


This  document  is  the  property  of  the  United  States  Government.  It  is  furnished  for  the  du¬ 
ration  of  the  contract  and  shall  be  returned  when  no  longer  required,  or  upon  recall  by  ASTIA 
to  the  following  address:  Armed  Services  Technical  Information  Agency, 
Document  Service  Center,  Knott  Building,  Dayton  2,  Ohio. 


i 

i 


v. 


NOTICE:  WHEN  GOVERNMENT  OR  OTHER  DRAWINGS,  SPECIFICATIONS  OR  OTHER  DATA 
Are  USED  FOR  ANY  PURPOSE  OTHER  THAN  IN  CONNECTION  WITH  A  DEFINITELY  RELATED 
GOVERNMENT  PROCUREMENT  OPERATION,  THE  U.  S.  GOVERNMENT  THEREBY  INCURS 
NO  RESPONSIBILITY,  NOR  ANY  OBLIGATION  WHATSOEVER;  AND  THE  FACT  THAT  THE 
GOVERNMENT  MAY  HAVE  FORMULATED,  FURNISHED,  OR  IN  ANY  WAY  SUPPLIED  THE 
SAID  DRAWINGS,  SPECIFICATIONS,  OR  OTHER  DATA  IS  NOT  TO  B  1  REGARDED  BY 
IMPLICATION  OR  OTHERWISE  AS  IN  ANY  MANNER  LICENSING  THE  HOLDER  OR  ANY  OTHER 
PERSON  OR  CORPORATION,  OR  CONVEYING  ANY  RIGHTS  OR  PERMISSION  TO  MANUFACTURE, 
USE  OR  SELL  ANY  PATENTED  INVENTION  THAT  MAY  IN  ANY  WAY  BE  RELATED  THERETO. 


UNCLASSIFIED 


il 


WPPlIW1 


% 


O 


Difference  Methods  on  a  Digital  Computer 
for  Laplacian  Boundary  Value 
and  Eigenvalue  Problems 

GEORGE  E.  FORSYTHE 


Reprinted  from 

COMMUNICATIONS  ON  PURE  AND  APPLIED  MATHEMATICS 
Volume  IX,  No.  3  August,  1956 

Copyright,  1958,  Lnterscienoe  Publishers,  Inc.  All  rights  reserved. 


COMMUNICATIONS  ON  PURE  AND 


APPLIED  MATHEMATICS,  VOL.  IX,  425-434  (1956) 


Difference  Methods  on  a  Digital  Computer 
for  Laplacian  Boundary  Value 
and  Eigenvalue  Problems* 

GEORGE  E.  FORSYTHE 

University  of  California ,  Los  Angeles 


1.  Finite  Difference  Approximations 

We  begin  with 

Definition  1.  A  simple  closed  curve  P:  x(s)-{-iy(s)  (0  ^  s  £2  st)  is  said 
to  be  piecewise  analytic  if  x(s)  and  y(s)  are  real  analytic  functions  of  the  arc 
length  s  of  F  in  each  of  a  finite  number  t  of  closed  intervals  0  =  s0  ^  s  ^  si  , 
sx  <:  s  <;  s2  ,  •  •  •  ,  s*-!  ^  s  ^  st . 

Let  the  open  region  R  of  the  x,  y-plane  have  a  boundary  C  which  is  the 
union  of  a  finite  number  of  piecewise  analytic  simple  closed  curves.  For  real- 
valued  functions  u  =  «(*,  y)  defined  on  R\jC,  denote  by  Au  the  Laplacian 
M**+M*»  .  In  numerical  computations  R,  C,  u,  and  Au  are  frequently  replaced 
by  corresponding  quantities  associated  with  a  finite  difference  network. 
These  may  be  introduced  as  follows  (see  [1]  for  a  similar  treatment): 

For  a  given  mesh  constant  h  >  0  let  a  net  consist  of  the  lines  x  =  fih, 
y  =  vh,  n,v  =  0,  ±1,  ±2-  ’  ’  *  •  The  points  (/ th ,  vh)  are  called  nodes.  The 
M  nodes  which  are  within  R  form  the  net  region  Rh  ,  assumed  to  be  connect¬ 
able  by  line  segments  of  the  net  within  R.  A  point  [jih ,  vh)  of  Rh  is  said  to  be 
a  regular  interior  point  if  each  of  the  four  neighboring  nodes  (jih±h,  vh±h) 
is  in  R  u  C.  All  other  points  of  Rh  are  called  irregular  interior  points.  Certain 
points  P  of  Rh  have  some  number  p  (1  2a  p  ^  4)  of  neighboring  nodes  Px  , 
.  .  .  t  pv  not  in  the  open  region  R.  For  any  such  point  P ,  let  Sk  denote  the 
closed  segment  of  the  net  joining  P  to  Pk  .  For  each  k  =  1,  •  •  •  ,  p,  the  point 
of  Sk  n  C  closest  to  P  will  be  called  a  boundary  point  of  the  net.  The  set  of 
boundary  points  will  be  denoted  by  Ch .  Some  points  of  Ch  may  be  nodes,  but 
each  irregular  interior  point  of  Rh  has  one  or  more  neighboring  points  in  Ch 
which  are  not  nodes.  (The  definition  of  Ch  given  in  [2],  although  intended  to 
be  the  same,  is  not  precise.) 


•The  research  leading  to  this  paper  was  sponsored  jointly  by  the  off“*  *e' 

search,  project  NR-044-144,  and  the  Office  of  Ordnance  Research,  project  TB*-0001(1210). 
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„  =  uS"5'"  “■  “ned  on  *UC,  is  a  real-valued  net-function 

For  each  point  °P  ^‘of  sTth  ®  *  d"""  d"finCd  “ 

denoted  by  (*+*,'  li  L  If  V I  L  4  "7,  ^  °‘  **uC *  be 

nf  *Hrv  7-  •  ^  E:y)’  (x'  y+h»)-  (*— ,  y),  and  (*,  «—/,_)  where  each 

of  the  A  s  is  positive  (£  ~  "east";  iV  ~  "north",  etc.)  Lei  777/7 

N  -ftJh,  etc.  For  regular  interior  points  bE  =  bN  = _ *1  Let 

Define  the  residual  rp  at  P  by 


(1)  >>  = 
where 


— —  — f— - f-  -  i _ MA7  % 


-4« 


p> 


B  -= 


46A 


+ 


4Ms‘ 


Then 

^  =  rP/h2. 

(The  definition  of  d<*>«  in  [2]  differs  by  the  factor  22?.)  In  case  P  is  a  reenlnr 
interior  point,  (1)  reduces  to  the  usual  formula  K 

^  rP  —  ue+uw-\-un-\-us~  Aup. 

When  there  are  irreguJar  interior  points>  the  tor  scJf.adjoint 


2.  Asymptotic  Inequalities  for  All  Eigenvalues 

For  the  operator  d,  the  eigenvalue  problem 


(4) 


d«+*M  =  0  in  R, 
u  =  0  on  C 

lSfna°nd  corrfuondin'f nit°  ,inearly  tadeP'"d«‘  eigenfunctions 

ana  corresponding  eigenvalues  *<*>,  numbered  so  that 

^  h  <  a(1>  <  A'*>  s;  a,3>  ^  . 

Similarly,  the  eigenvalue  problem 

(6)  ^am+^a«  =  0  on  1?,, 

«  =  0  on  Ch 

blues' Ar.Xr  °' ,incarly  indeIX'nden,  ^“-‘‘ons  „■»  and  eigen- 

0  <  '-I11  <  ^  2<8>  ^  g  ;.w. 

(Recall  that  JW  is  the  number  of  points  of  Rh  .) 
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For  rectangles  whose  sides  are  lines  of  the  net,  it  is  well  known  that 
(7)  <  X-’.  k  =  1,  •  •  • ,  M; 

see  also  below.  Since  for  general  regions  R  close  lower  bounds  for  the  A(fc) 
are  not  easily  found,  the  author  has  been  investigating  [3,  5a,  2,  4]  whether 
(7)  might  also  hold  for  other  regions. 

Suppose  one  deals  with  a  polygonal  region  R  which  is  a  finite  union  of 
squares  and  half-squares  of  a  net  of  mesh  constant  2 h.  Then  R^  has  no  ir¬ 
regular  interior  points,  and  neither  does  Rh  ,  obtained  by  refining  the  net. 
At  any  point  ( x ,  y)  of  Rn  , 

(®)  ^2 =  \bfiL,ylU, 

where  Lhu  =  h~*{uxx-x-+u„~}.  Here 

uxx5x  =  u{x—2h,  y)—4u(x—h,  y)+6u(x,  y)—4u(x+h,  y)+u(x+2h,  y), 

and  uvvii  is  defined  analogously.  Following  the  idea  [3]  of  putting  the  fun¬ 
damental  eigenfunction  uh  =  ujf1  of  the  problem  for  Rh  into  the  minimization 
problem  for  RM  ,  one  finds  that 

<»>  Ssiy-ti’EVA)  (2  «»-*. 

(2*)  (2A) 

where  denotes  a  sum  over  the  points  of  R ^  .  Hence  whenever 
(2*) 

(10)  >  o 

(2*) 

it  follows  that 

(11)  <A«o. 

If  (11)  should  hold  also  for  h  replaced  by  A/2,  A/4,  A/8,  •  •  • ,  it  would  then 

follow  that  A^1’  <  A(1),  since  the  methods  of  [5]  show  that  Aj1*  ->  A(1>,  as 

A  ->  0. 

Recent  computation  with  SWAC  (see  Section  4  below)  for  the  L- 
shaped  region  of  [3]  shows  that  (11)  fails  for  A  =  1/16.  Hence  (10)  is  not 
true  for  all  polygonal  regions  of  the  type  considered  here.  I  conjecture, 
however,  that  (10)  holds  for  any  convex  R  which  is  the  union  of  squares  and 
half-squares  of  the  net. 

Lacking  a  proof  of  this  conjecture,  the  author  proved  [3]  an  asymptotic 
form  of  (7)  for  k  =  1,  and  in  [5a]  and  [2]  generalized  it  to  include  certain 
regions  with  curved  boundaries  (the  case  k  =  1  of  Theorem  2  below).  Thus, 
under  the  hypotheses  of  Theorem  2  below,  the  desired  inequality  (7)  has  al¬ 
ready  been  proved  for  the  first  eigenvalue,  though  only  for  sufficiently  small 
A.  Less  delicate  inequalities  of  type  Aj1J  <  A*1’,  valid  for  all  A,  have  since 
been  found  by  Weinberger  [6,  7]  and  Hersch  [8]  without  such  stringent 
hypotheses  on  R,  by  allowing  Rh  to  correspond  to  a  region  a  little  larger  than 
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R.  Assuming  *(s),  y(s)  defining  r  to  be  seven  times  continuously  differen¬ 
tiable,  Saul’ev  [8a]  has  recently  found  inequalities  of  the  type  l^’— yt(i)|  = 
0(h2~‘),  as  h  ->  0. 

During  the  1955  Berkeley  Conference  for  Partial  Differential  Equations, 
at  which  this  paper  was  presented,  the  author  completed  an  extension  of  the 
results  of  [2]  to  all  eigenvalues  of  A  and  Ah  .  This  has  been  abstracted  in  [4], 
and  will  now  be  summarized.  Let  r  be  ihe  angle  between  the  *-axis  and  the 
positive  direction  of  the  tangent  to  C.  For  any  function  v,  let  vn  denote  the 
inner  normal  derivative  of  v  on  C.  For  any  two  functions  v,  w,  define  an 
inner  product  ( v ,  w)  by  the  formula 

12(v,  w)  =  jj(vxxwxx-\-vvvwl/v)  dx  dyJrjvnwn  sin2  2r  dx. 

R  c 

Let  the  infinite  symmetric  matrix  D  have  elements  d(i  =  («<*>,  uu>),  i,  j  — 
1,  2,  •  •  •  ,  where  u{k)  is  an  eigenfunction  of  (4)  belonging  to  A<*>,  normalized 

so  that  //.  (ulk))2dx  dy  =  1.  Corresponding  to  an  m-fold  multiple  eigenvalue 
of  A,  we  are  to  take  m  linearly  independent  functions  u{k). 

Define  a  sequence  of  real  numbers  y<*>,  k  =  1,  2,  •  •  •,  as  follows.  If 
A1*'  is  a  simple  eigenvalue  of  (4),  take  y(l>  =  dkk  .  If  A(fc>  =  =  .•■== 

^<*+m-i)  js  an  ^.foid  multiple  eigenvalue  of  (4),  then  define  y<fc>  5>  y(*+i>  ^ 

.  .  .  ^  yik+m-i)  to  m  (rea^  eigenvalues  of  (d(j),  i,  j  =  k,  •  •  •  ,  k-\-m— 1, 

the  corresponding  w-rowed  principal  minor  of  D. 

definition  2.  A  region  R  bounded  by  a  finite  number  of  piecewise 
analytic  curves  is  said  to  have  strictly  convex  corners  if,  at  any  point  x(s,)  + 
iy(Sj)  where  distinct  analytic  curves  meet  ( see  Definition  1),  the  interior  angle  is 
less  than  ji  ( equality  not  permitted). 

theorem  1.  Let  R,  an  open  region  with  strictly  convex  corners,  be  bounded 
by  C,  the  union  of  a  finite  number  of  piecewise  analytic,  simple  closed  curves. 
Define  h,  A(k>  ,  ffl,  y{k)  as  above.  Then  for  k  =  1,  2,  •  •  •  one  has  —  oo  <  < 

oo  and,  as  h  ->  0, 

(12)  ^  A‘*>-y<*>A2+o(/i2). 

A  proof  of  Theorem  1  will  appear  in  a  later  paper.  If  R  is  also  convex, 
then  C  is  a  simple  closed  curve  with  dx  ^  0,  and  D  is  therefore  a  Gramian 
matrix.  Hence  all  y{k)  >  0,  and  we  have  the  following  corollary  to 
Theorem  1: 

theorem  2.  Let  R,  a  convex  open  region  with  strictly  convex  corners,  be 
bounded  by  C,  a  piecewise  analytic,  simple  closed  curve.  Then,  with  the  above 
notation,  0  <  y<k>  <  oo,  and  for  each  k  =  1,  2,  •  •  •  there  exists  h0(k)  >  0  such 
that  <  A<*>  if  h  ^  h0{k). 
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The  case  k  —  1  of  Theorems  1  and  2  is  the  result  proved  in  [2].  (In 
Theorem  1  of  [2]  the  assumption  that  R  be  simply  connected  was  never 
used.  Note:  The  proof  of  Lemma  1  in  [2]  is  incorrect.) 

As  an  asymptotic  result  Theorem  1  is  the  best  possible,  in  the  following 
sense.  For  a  rectangle  R  of  sides  nfp,  nfq,  the  relation  (12)  is  actually  an 
equality  for  all  k,  up  to  the  term  o(h2).  For,  corresponding  to  the  eigenfunc¬ 
tion  «<"■•>  =  sin  mpx  sin  nqy  of  both  (4)  and  (6)  (it  is  convenient  here  to 
use  a  double  superscript  m,  n  instead  of  the  single  superscript  k),  one  has 

=  {mp)*+{nq)\ 


and,  as  h 


(13) 


0  over  values  of  h  for  which  R  is  a  union  of  squares  of  the  net, 
**n,n’  =  sin2  (ntph/2)  +  sin2  (nqh/2) 

( W  ” 


=  yP*»,n)  _ 


{mp)*-\-  ( nq )* 
12 


h2+o(h2) 


■  u{m-n>)h2/12+o(h2). 

Since  (u(  \  u^-  ()  =  0  imless  \tn—(i\-\-\n—v\  =  0,  one  has  — 

(m  *»-»),  M(”>.«))/i2(  even  for  multiple  eigenvalues.  Hence  (13)  shows  that 
(12)  is  an  equality  for  the  rectangle  R,  up  to  the  term  o(h2). 


3.  Solving  the  Difference  Equations  on  a  Machine 


Let  us  first  consider  solving  the  Dirichlet  difference  problem 
(14a)  A*-h)u  =  0  on  Rh, 

(14b)  «  =  prescribed  function  on  C. 

n 


on  an  automatic  digital  computer.  Assume  there 
points  of  Rh  ,  so  that  (14a)  is  equivalent  to 


are  no  irregular  interior 


(14c)  rP  ~  ue+un+Uw+Us—*Up  =  0  on  Rh  . 

To  solve  the  system  of  M  equations  (14c)  in  the  unknown  values  of  u  at  the 
US,°f  **'  '\erltive  methods  are  usually  preferable  to  elimination 
^  n  ni  °f  th6ir  simPlicity-  The  "«««**  overrelaxation  process 

of  [10]  and  [9]  has  proved  very  successful.  In  it  the  points  of  Rh  are  re¬ 
peatedly  scanned  in  a  cyclic  order.  In  its  turn,  at  each  point  P  of  R„  the 
value  up  is  immediately  replaced  by  a  new  value  u'p , 

up  =  up~hipPP  , 


where  (1  ^  /S  <  2)  is  a  parameter  to  be  discussed  below, 
value  of  tip  ,  one  goes  on  to  the  next  point  of  Rh  ,  and  so  on. 


Using  the  new 
For£  =  1  this 
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» 


is  the  Liebmann  process  [11],  in  which  u'P  =  \(uE+uN+uw+us).  Any 
reasonable  initial  va’ues  uf]  of  the  function  u  over  Rh  will  suffice,  and  one 
may  take  all  uf]  =  0. 

Let  ujik)  be  the  result  of  k  sweeps  through  Rh  with  successive  overrelaxa¬ 
tion.  Let  «(®>  be  the  exact  solution  of  (14);  define  the  error  ef*  =  M(®>. 

Since  the  solution  process  is  linear,  one  knows  that 

4*  =  Hkef>  , 

where  H„  is  a  certain  linear  transformation  of  RK  into  itself.  The  conver¬ 
gence  of  ej,k)  to  0,  as  k  -*»  oo,  depends  on  the  maximum  modulus  of  the  eigen¬ 
values  of  Hp  .  For  the  Dirichlet  problem  it  is  known  [9]  that  |Al(f//})| 

<  1  (all  i,  0  <  fi  <  2),  so  that  successive  overrelaxation  (1  rg  p  <  2)  always 
converges.  The  ultimate  speed  of  convergence  depends  on  max,  \Xt(Hp)\. 
Now,  for  any  single  step  (15),  to  choose  p  =  1  is  the  best  tactic,  in  the  sense 
that  the  negative  definite  error  function  ^RcAwe  is  brought  as  close  to  the 
desired  value  0  as  is  possible  in  the  one  step,  by  choosing  p  =  1  in  (15). 
However,  hand  computers  have  long  observed  [12]  that  certain  choices 
P  >  1  apparently  represent  a  better  strategy  than  P  =  1,  in  the  sense  that 
e(pk)  -*■  0  faster  asymptotically,  as  k  -►  oo.  (Related  ideas  of  strategy  and 
tactics  for  other  linear  processes  will  be  found  in  [13].)  In  [10]  and  [9]  the 
matter  of  a  best  strategy  is  settled  by  showing  how  to  pick  the  unique  po 
for  which  max,  |A,(//j)|  is  minimized.  Indeed, 


(16) 

where 


^0  i  i  / :  ■  > 

1  + V  1  —  a 


(17) 


a  —  max,  |A,(^j)|. 


The  number  a  is  easily  approximated  by  setting  p  =  1  and  estimating  the 
asymptotic  value  of  ||^*+1)||/||ci,:,||,  as  k-r  oo,  where  |(-|j  denotes  some  con¬ 
venient  norm  function.  When  p  =  p0  ,  all  the  |A,(/^)|,  i  =  1,  •  •  •  ,  Af,  turn 
out  to  be  equal,  so  that  the  manner  of  approach  of  ||^w||  to  0  is  utterly  un¬ 
predictable. 

When  R  has  a  curvilinear  boundary  C,  the  theory  of  [9]  still  holds,  and 
the  successive  overrelaxation  method  is  just  as  attractive.  One  simply  re¬ 
places  the  formula  (3)  for  rp  at  irregular  interior  points  P  of  Rh  by  the  for¬ 
mula  (1). 

For  the  eigenvalue  problem  (6)  successive  overrelaxation  is  readily 
adaptable  to  getting  Aj1*  and  its  associated  eigenvector  ujf*  .  One  uses  a  trial 
value  A  for  Aj;11,  and  solves  the  system 


(4*+A/)m  =  0  on  Rh, 
u  —  0  on  C, 


L 


(18) 


COMPUTING  LAPLACIAN  PROBLEMS 


431 


by  successive  overrelaxation  starting  with  some  nonzero  vector  «i°!  •  For  the 
residual,  assuming  that  Rh  has  only  regular  interior  points  (so  that  Ah  = 
one  uses 


rP  —  W£  +  WJV  +  MH/  +  %  —  4:Up-{-AUp 

in  (15).  If  X  were  exactly  equal  to  X™  ,  the  iterates  uf  would  stay  roughly 
m  the  affine  subspace  S  (of  Af-dimensional  space)  which  is  orthogonal  to 
uh  and  passes  through  uf]  .  Apart  from  an  exceptional  set  of  «i°>  of  zero 
measure,  one  would  always  have  «<*>  ->  ««»> ,  as  k  ->  oo,  where  «  is  some 
nonzero  number  depending  on  u<°>  .  If  X  is  near  X™  but  not  equal  to  it,  one 
finds  that  u'g  moves  asymptotically  into  the  line  xu™  while  converging  to 
0  or  oo.  An  improved  value  of  X  is  then  obtainable  from  the  Rayleitrh  quo¬ 
tient  eh «’),  where  1 


One  then  uses  the  better  X  in  (19),  and  alternates  a  few  steps  of  succes¬ 
sive  overrelaxation  with  the  Rayleigh  quotient  calculation,  until  a  satis¬ 
factory  a is  obtained. 

/^or  Siting  Xlk>  and  «<*>  (k  >  1),  the  idea  of  solving  (18)  with  X  near 
Xh  is  undoubtedly  still  a  good  one,  but  successive  overrelaxation  no  longer 
converges,  because  Ah+X^I  is  not  semidefinite  for  k  >  1.  A  simple  iteration 
which  can  solve  indefinite  systems  must  be  substituted— perhaps  that  of 
Kaczmarz  [14], 


4.  SWAC  Codes 

SWAC  is  an  electronic  digital  computer  built  by  the  National  Bureau  of 
Standards,  but  now  owned  by  the  Department  of  Defense  and  loaned  to  the 
University  of  California,  Los  Angeles,  for  research  in  numerical  analysis. 
It  is  a  three  or  four  address  machine  with  a  Williams  tube  memory  of  256 
words;  each  word  consists  of  36  binary  digits  and  a  sign.  Add  time  is 
.000064  second;  multiplication  time  is  slightly  less  than  six  add  times.  An 
auxiliary  magnetic  drum  memory  holds  4096  words  in  128  channels  of  32 
words  each  (increased  to  8192  words  in  late  1955).  One  can  transfer  a  block 
of  up  to  32  words  between  the  drum  and  the  Williams  tube  memory  in 
.017  second.  For  substantial  problems,  computing  time  is  measured  in 
hours,  while  coding  time  is  measured  in  weeks. 

In  connection  with  the  author  s  graduate  seminar  on  numerical  methods 
for  elliptic  partial  differential  equations,  held  in  the  spring  of  1955,  a  few 
general  codes  were  prepared  following  the  methods  of  Section  3  above. 

In  all  these  codes  the  net  region  Rh  is  represented  by  a  subset  of  the 
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32  by  128  rectangular  array  of  cells  on  SWAC’s  magnetic  drum.  Parameters 
indicate  the  smallest  rectangle  in  which  Rh  is  contained.  In  all  the  codes  the 
data  in  successive  drum  channels  are  brought  to  the  Williams  tube  memory 
for  the  relaxation  process,  and  the  improved  values  are  sent  back  when  they 
are  no  longer  needed. 

The  author  s  basic  Code  60  solves  Dirichlet’s  problem  for  the  operator 
A(h)  ==  Ah  for  a  region  with  no  irregular  interior  points.  Here  Ch  consists 
solely  of  nodes  of  the  net,  and  can  thus  also  be  represented  by  cells  of  the 
drum.  Points  of  Rh  have  a  0  in  the  most  significant  binary  digit.  Points  not 
in  Rh  are  tagged  with  a  1  in  the  same  digit,  and  can  be  an  arbitrary  subset  of 
the  32  by  128  rectangular  array.  Code  60  assumes  that  the  values  of  u  have 
been  placed  in  the  drum  cells  representing  Ch  ,  in  digits  less  significant  than 
the  first  (which  is  assumed  to  be  tagged  with  a  1,  as  stated  above).  It  per¬ 
forms  successive  overrelaxation  at  all  points  of  Rh,  leaving  all  other  values  of 
u  unchanged.  The  overrelaxation  factor  £  defined  in  Section  3  can  be  modi¬ 
fied  at  will  by  the  SWAC  operator.  (It  is  fun  for  the  operator  to  participate 
in  the  solution.  Having  the  operator  manage  /?,  while  SWAC  manages  all 
other  numbers,  seems  a  satisfactory  division  of  labor!)  The  error  measure 
used,  denoted  by  ||^>||,  is  the  first  power  norm  of  r<*>;  it  is  displayed  at  the 
end  of  each  sweep  through  Rh  .  It  is  easily  accumulated  during  the  algorithm 
(15),  since 

Ill'll  -  X  l-W 

The  ratio  ||4*,||/||*jt*“1,||  is  also  displayed,  and  is  useful,  for  example, 
in  approximating  the  a  of  (17). 

When  Rh  is  a  rectangle  of  30  by  68  interior  points,  a  relaxation  sweep 
through  the  2040  points  takes  SWAC  8.5  seconds.  When  /?  =  1,  to  reduce 
by  a  factor  of  10~8  would  take  about  2300  sweeps — over  five  hours  on 
SWAC.  Using  the  optimal  value  /5  ==  /50  =  1.858-f-,  only  90  sweeps  would  be 
needed  about  13  minutes  on  SWAC.  Including  the  time  necessary  to 
estimate  that  1.875  is  a  pretty  good  value  for  /5,  and  running  with  this  /?,  a 
reduction  of  10~8in  e{pk)  is  found  in  practice  to  take  20  or  30  minutes  (see  [15]). 

Code  61  is  a  modification  of  Code  60  in  which  boundary  conditions  cor¬ 
responding  to  du/dn  =  un  =  0  are  permitted  at  any  point  of  Ch  .  Suppose, 
for  example,  that  it  is  desired  to  have  dujdx  =  0  for  x  =  0,  where  R  is  a 
region  with  x  >  0.  Then  values  of  u(x,  y)  for  the  lines  x  =  —h,  x  =  0, 
x  —  h,  x  —  2 h,  •  •  •  are  all  stored  on  the  drum  for  appropriate  values  of  y. 
When  the  code  comes  to  «( — h,  y)  for  any  y,  the  value  is  simply  replaced  by 
u(h>  y)-  The  values  u( 0,  y),  u(h,  y),  u(2h,  y),  •  •  •  are  all  relaxed  according  to 
(15).  Again  the  relaxation  routine  is  notified  of  the  different  types  of  interior 
and  boundary  points  by  suitable  identification  digits. 
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Mr.  Kenneth  Ralston  of  U.C.L.A.  has  prepared  a  SWAC  code  analogous 
to  Code  60,  dealing  with  curvilinear  boundaries.  Note  in  this  case  that  it  is 
no  longer  possible  to  map  Ch  onto  the  drum  in  a  simple  manner.  For  example, 
the  point  of  Ch  "north"  of  the  point  (fih,  vh)  of  Rh  may  not  be  the  same  as  the 
point  of  Ch  "west”  of  ( ,  vh-\-h).  Identification  digits  distinguish: 
(a)  points  not  in  Rh  ,  to  be  left  alone,  (b)  regular  interior  points,  to  be  treated 
by  (15),  with  rP  as  in  (3),  and  (c)  irregular  interior  points,  to  be  treated  by 
(15),  with  rp  as  in  (1). 

To  deal  with  irregular  interior  points  Ralston  stores  the  coefficients  of 
11  e  >  uw  >  un  >  «s  *n  formula  (1)  in  four  cells  of  certain  special  channels  of  the 
drum.  Since  one  meets  the  irregular  interior  points  of  Rk  in  a  fixed  order,  all 
these  special  coefficients  can  be  stored  once  and  for  all  in  the  proper  sequence. 
(Such  of  the  uE  ,  uw ,  uN  ,  us  as  arc  fixed  boundary  values  are  actually 
multiplied  into  the  coefficients  stored  in  the  special  channels.)  Ralston  can 
also  use  the  special  channels  to  simulate  a  general  mixed  boundary  condition 
<*un+bu+c  =  0  on  a  straight  boundary  segment  parallel  to  the  x  or  y  axis. 

The  author’s  Code  70  performs  successive  overrelaxation  for  the  eigen¬ 
value  problem  (6)  over  regions  with  no  irregular  interior  points,  following 
the  method  described  in  Section  3.  The  quantities  /?  and  4— A  arc  both 
parameters  of  the  code.  Auxiliary  to  Code  70  are  three  other  codes:  Code  80 
doubles  the  values  of  u  over  Rk  ,  while  Code  82  halves  these  values  of  u. 
Code  81  forms  the  Rayleigh  quotient  (20). 

With  Codes  70—82  the  author  found  the  fundamental  eigenvalue  for 
the  /.-shaped  region  of  [3]  for  various  h,  and  found  the  counter-example 
mentioned  after  (11).  Code  61  has  been  adapted  for  some  soil  engineering 
calculations  by  Mr.  D.  Isherwood  and  Mrs.  Louise  Straus  of  U.C.L.A.  Code 
60  has  been  modified  by  Miss  Rita  Powers  at  U.C.L.A.  to  study  the  influence 
of  the  order  in  which  the  points  of  Rh  arc  relaxed  on  the  rate  of  convergence 
of  successive  overrelaxation. 

To  make  Ralston’s  code  more  available  for  problems,  a  code  is  badly 
needed  to  tag  the  drum  cells  with  the  identification  digits,  to  precompute  the 
coefficients  (1),  and  to  store  the  latter  on  the  special  channels.  Ralston  and 
the  author  have  outlined  such  a  code,  assuming  C  to  be  a  piecewise  quadratic 
function  of  x,  y.  Another  needed  code  is  one  to  get  A[w,  for  k  >  1. 

Using  the  ORACLE  computer  at  Oak  Ridge  National  Laboratories, 
W.  C.  Sangren  [unpublished]  has  also  computed  AJ11  for  the  L-shaped  region 
of  [3]  and  for  other  simple  regions  with  rectilinear  boundaries. 
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